﻿using System;
using System.Collections.Generic;
using System.Text;
using System.Collections;
using VerticalProcesses;
namespace VerticalProcesses
{
    /// <summary>
    /// Basic snowpack unit
    /// </summary>
    /// 
    public abstract class SnowpackType
    {
        // state and flux variables
        protected double Mf = 0.0;                     // potential snow melt rate, 
        protected double RainSnow = 0.0;               // rain on snow
        protected double SweSurf = 0.0;                // Water equivalent of snow surface layer, mm
        protected double SwePack = 0.0;                // Water equivalent of snow pack layer, mm
        protected double SlwcSurf = 0.0;               // Liquid water content of snow surface layer, mm
        protected double SlwcPack = 0.0;               // Liquid water content of snow pack layer, mm


        // parameters
        protected double MAXLIQ = 0.05;                // maximum liquid water fraction of snow, dimensionless
        protected double MAXSWESURF = 100;             // maximum swe of surface layer of snowpack

        #region public properties

        public double Swe {get;set;}
        public double SnowOutflow { get; set; }
        public double TempSurf {get;set;}
        public double TempPack {get;set;}
        public double CcSurf {get;set;}    
        public double CcPack {get;set;}

        #endregion

        /// <summary>
        /// updates snowpack for a given timestep
        /// </summary>
        public abstract void Update(double snowpackEnergyInput, double thrfRain, double thrfSnow, double gruTa);

    }

    # region Implementions of different Snowpacks
    /// <summary>
    /// Add your snowpack implemention here
    /// </summary>
    public class SnowpackTwoLayer : SnowpackType
    {


        /// <summary>
        /// Snowpack with two layers
        /// surface layer exchanges engergy with athmosphere
        /// pack layer gets water from surface layer
        /// </summary>
        /// <param name="snowpackEnergyInput"></param>
        /// <param name="thrfRain: rain throughfall rate [mm/h]"></param>
        /// <param name="thrfSnow: snow throughfall rate [mm/h]"></param>
        /// <param name="gruTa: air temperature [dC]"></param>

        public override void Update(double snowpackEnergyInput, double thrfRain, double thrfSnow, double gruTa)
        {

            //initialize flux varaibles
            SnowOutflow = 0.0;
            Mf = 0.0;
            //IN and OUT
            //cc                      cold content of snowpack (positive), MJ/m2

            //OUPUT
            //snowMelt               melt rate from snowpack, mm/h
            //_rainSnow               rain added to snowpack, mm/h


            //CONSTANTS (SEE CONSTANTS)
            //CvIce - volumetric heat capacity of ice, MJ m-2 mm-1 K-1
            //LF  heat of fusion of water, MJ m-2 mm-1

            //INTERNAL
            double avLiq;               //available liquid storage of snow
            double meltSurf = 0.0;      //melt rate from snow surface layer, mm/h  
            double posMf;               // positive meltfactor
            double ccSnowFall;

            //reconstruct two layer snowpack
            //------------------------------------------------------------------------
            if (Swe > MAXSWESURF)
            {
                SweSurf = MAXSWESURF;
            }
            else
            {
                SweSurf = Swe;
            }
            SwePack = Swe - SweSurf;

            //add snow throughfall to snowPack and compute its cold content
            //this also initializes the cold content with the first snow fall
            //------------------------------------------------------------------------

            CcSurf = Constants.CvIce * SweSurf * -TempSurf;
            CcPack = Constants.CvIce * SwePack * -TempPack;
            ccSnowFall = Constants.CvIce * thrfSnow * Math.Max(-gruTa, 0.0);

            //distribute fresh snowfall
            //------------------------------------------------------------------------

            if (thrfSnow < (MAXSWESURF - SweSurf))// suface layer is below max
            {
                SweSurf += thrfSnow;
                CcSurf += ccSnowFall;
            }


            else
            {//snow surface layer is at max
                double deltaswePack = SweSurf + thrfSnow - MAXSWESURF;
                double deltaccPack = deltaswePack / SweSurf * CcSurf;//give pack the proportional coldconent of surface layer
                SwePack += deltaswePack;
                SweSurf = MAXSWESURF;
                CcPack += deltaccPack;
                CcSurf += ccSnowFall - deltaccPack;
            }

            //freeze liquid water in surface layer
            //cant have cold content and liquid water-one needs to be 0
            //------------------------------------------------------------------------

            if (CcSurf > 0.0 & SlwcSurf > 0.0)
            {
                if (CcSurf > SlwcSurf * Constants.LF)
                {
                    CcSurf -= SlwcSurf * Constants.LF;        //freeze all liquid water
                    SlwcSurf = 0.0;
                }
                else
                {
                    SlwcSurf -= CcSurf / Constants.LF;      //freeze only part of water
                    CcSurf = 0.0;
                }
            }


            //snowpack cooling or warming
            //engergy exchange occurs only between snow surface layer and atmosphere
            //snow surface layer and snow pack layer only exchange liquid water
            //------------------------------------------------------------------------

            //Convert melt energy [MJ/m2/timestep] to melt rate [mm/h]
            //negative _Mf is cold content in [mm/h] (later converted back to [MJ/m2])
            //------------------------------------------------------------------------

            if (SweSurf > 0.0)
            {
                Mf = snowpackEnergyInput / Constants.LF;

                //snowpack cooling or warming
                //------------------------------------------------------------------------
                if (Mf <= 0.0)
                {      //   snow cooling
                    posMf = -Mf;
                    if (posMf < SlwcSurf)//refreeze part of liquid water
                    {
                        SlwcSurf = SlwcSurf - posMf;        // note that _Mf is negative (decrease _slwcSurf)
                        CcSurf = 0.0;
                    }
                    else// refreeze all liquid water and increase cold content
                    {
                        posMf = posMf - SlwcSurf;
                        SlwcSurf = 0.0;
                        CcSurf = CcSurf + posMf * Constants.LF;     //convert snowpackEnergyInput back to MJ/m2

                        if (SweSurf < 10.0) //dont cool below gruTa when Swe is below 10 mm (need this for stable run, otherwise CC and Tsnow start to fluctuate...)
                        {
                            CcSurf = Math.Min(CcSurf, -gruTa * SweSurf * Constants.CvIce);
                        }

                    }
                }
                else
                {     // snow warming (_Mf > 0)
                    if (Mf * Constants.LF < CcSurf | gruTa < 0.0)
                    {
                        //not enough energy to satisfy cold content at surface layer and to start melt
                        //reduce cold content and step out
                        if (gruTa < 0.0)
                        {
                            CcSurf = CcSurf - Mf * Constants.LF;
                            //dont warm above gruTa                           
                            CcSurf = Math.Max(CcSurf, -gruTa * SweSurf * Constants.CvIce);
                        }

                        else
                        {
                            CcSurf = CcSurf - Mf * Constants.LF;
                        }

                        SlwcSurf = 0.0;      //should be 0 already...
                    }

                    else
                    {   //eliminate CC (if there is any)
                        Mf = Mf - CcSurf / Constants.LF;
                        CcSurf = 0.0;

                        if (Mf <= MAXLIQ * SweSurf - SlwcSurf)   //all snow melt goes into liquid water storage
                        {     //still liquid water storage left
                            SlwcSurf = SlwcSurf + Mf;
                        }
                        else
                        {     //liquid water storage full, snow melt...
                            Mf = Mf - (MAXLIQ * SweSurf - SlwcSurf);      // reduce the available melt by the part the goes to liquid water storage

                            if (SweSurf * (1.0 - MAXLIQ) > Mf)  // not enough energy to melt entire surface layer
                            {
                                meltSurf += Mf / (1.0 - MAXLIQ);      //snow plus liquid in it
                                SweSurf -= Mf / (1.0 - MAXLIQ);
                                SlwcSurf = MAXLIQ * SweSurf;
                            }
                            else
                            {
                                meltSurf += SweSurf;      //all snow melts plus plus the rainfall
                                SweSurf = 0.0;
                                SlwcSurf = 0.0;
                            }
                        }

                    }
                }
                //rain on snow surface layer
                //note that the rain energy input is handled in snowpackengery input
                //snowpackenergy input was handled before
                //as a result there is no melt here but only increase of swe
                //------------------------------------------------------------------------
                if (thrfRain == 0.0 | SweSurf == 0.0)
                {
                    RainSnow = 0.0;
                    meltSurf+=thrfRain;
                }
                else
                {
                    if (CcSurf > 0.0)
                    {   //CC refreezes rain
                        //engery input of rain is already added in snowpackEnergyInput
                        if (CcSurf > thrfRain * Constants.LF)
                        {     //refreeze all rain
                            CcSurf = CcSurf - thrfRain * Constants.LF; //derease CC
                            RainSnow = thrfRain;
                            SweSurf = SweSurf + thrfRain;
                        }
                        else
                        {     //refreeze only part of rain and eliminate CC
                            RainSnow = CcSurf / Constants.LF;
                            SweSurf = SweSurf + RainSnow;
                            CcSurf = 0.0;
                            //remaining rain
                            thrfRain = thrfRain - RainSnow;      //reduce rain by part that refroze
                            
                            // now check if there is still some liquid storage left
                            if (thrfRain < MAXLIQ * SweSurf / (1.0 - MAXLIQ))//still some liquid storage left
                            {     //remaing rain becomes _slwcSurf
                                SlwcSurf = thrfRain;
                                RainSnow += thrfRain;       //add liquid part back
                                SweSurf += thrfRain;
                            }       //add liquid part
                            else // no liquid storage left
                            {
                                SlwcSurf = MAXLIQ * SweSurf / (1.0 - MAXLIQ);
                                RainSnow += SlwcSurf;
                                SweSurf += SlwcSurf;
                                meltSurf = meltSurf + thrfRain - SlwcSurf;
                            }
                        }
                    }
                    else
                    {     //CC =  0
                        if (SlwcSurf >= MAXLIQ * SweSurf)
                        {     //_MaxLiquid already reached (all rain on snow goes into meltsurf)
                            RainSnow = 0;
                            meltSurf += thrfRain; 
                        }
                        else
                        {
                            avLiq = MAXLIQ * SweSurf - SlwcSurf;
                            if (thrfRain < avLiq)
                            {     // all rain into snow (this._slwcSurf)
                                SlwcSurf += thrfRain;
                                RainSnow = thrfRain;
                                SweSurf += thrfRain;
                            }
                            else
                            {     //_MaxLiquid reached)
                                RainSnow = avLiq / (1 - MAXLIQ);
                                SweSurf += RainSnow;
                                SlwcSurf = MAXLIQ * SweSurf;
                                thrfRain -= RainSnow;
                                meltSurf += thrfRain;
                            }
                        }
                    }
                }
            }
            //snow pack layer
            //------------------------------------------------------------------------
            //model runs as a single layer melt model when in absence of a snow pack layer
            //the second layer (snowpack) gets its cold content from the suface layer
            if (SwePack > 0.0)
            //add the surface layer outflow to pack liquid water 
            {
                SlwcPack += meltSurf;
                //SwePack += meltSurf;

                //re-freeze liquid water in snowpack layer
                //------------------------------------------------------------------------
                if (CcPack > SlwcPack * Constants.LF) //snow pack layer cold content not fully satisfied
                {
                    CcPack -= SlwcPack * Constants.LF;        //freeze all water
                    SlwcPack = 0.0;
                    SwePack += meltSurf; //all water gets added to pack layer
                }
                else //eliminate cold content
                {
                    SlwcPack -= CcPack / Constants.LF;      //freeze only part of water
                    SwePack += meltSurf; //part of the water stays in liquid storage but SWE is the total 
                    CcPack = 0.0;
                    
                    //update snow pack layer liquid water content
                    //------------------------------------------------------------------------
                    if (SlwcPack > SwePack * MAXLIQ)//liquid water storage full, start outflow               

                    //note that if swePack == 0,i.e., when model switches to one layer<
                    //swePack * MaxLiquid = 0 and snowOutflow = meltSurf
                    {
                        SnowOutflow = SlwcPack - SwePack * MAXLIQ;
                        SlwcPack = SwePack * MAXLIQ;
                        SwePack -= SnowOutflow;
                        SwePack = Math.Max(0.0, SwePack);
                    }
                }


            }
            else // not pack layer
            {
                SnowOutflow = meltSurf;
            }
            // clean up
            Swe = SweSurf + SwePack;
            TempSurf = -CcSurf / (Constants.CvIce * Math.Max(1, SweSurf));//this assumes an isothermal snow surface layer
            TempPack = -CcPack / (Constants.CvIce * Math.Max(1, SwePack));//this assumes an isothermal snow surface layer
        }
        public SnowpackTwoLayer(double MAXLIQ, double MAXSWESURF, double IniSwe)
        {
            this.MAXLIQ = MAXLIQ;
            this.MAXSWESURF = MAXSWESURF;
            this.Swe = IniSwe;
        }

    }
}
#endregion